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In this paper, we study triply periodic surfaces with minimal surface area under 
a constraint in the volume fraction of the regions (phases) that the surface separates. 
Using a variational level set method formulation, we present a theoretical charac- 
terization of and a numerical algorithm for computing these surfaces. We use our 
theoretical and computational formulation to study the optimality of the Schwartz 
P, Schwartz D, and Schoen G surfaces when the volume fractions of the two phases 
are equal and explore the properties of optimal structures when the volume fractions 
of the two phases not equal. Due to the computational cost of the fully, three- 
dimensional shape optimization problem, we implement our numerical simulations 
using a parallel level set method software package. 



1. Introduction 

Triply periodic minimal surfaces [111121 El are objects of great interest to physical 
scientists, biologists, and mathematicians because they naturally arise in a variety 
of systems, including block copolymers [E], nanocomposites [El, micellar materials [ 
ini, and lipid- water systems and certain cell membranes [ 13111013 ■ There are two 
key feature of these surfaces: (1) the mean curvature^ is zero everywhere on the 

"'^The mean curvature, H(x.), at a point x of a surface is the average of the two principal normal 
curvatures, Ki(x) and K2(x): H{x) = ^ +K2(x)). i?(x) is also conveniently represented in 
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surface and (2) they are periodic in all three coordinate directions (i.e. they extend 
infinitely in all directions and possess the symmetry of one of the crystallographic 
space groups). An important subclass of triply periodic minimal surfaces are those 
that partition space into two disjoint but intertwining regions that are simultaneously 
continuous [ 10^. Examples of such surfaces include the Schwartz primitive (P), the 
Schwartz diamond (D), and the Schoen gyroid (G) surfaces (see Figure HJ. 




Figure 1. Unit cell for three common minimal surfaces: Schwartz P surface (left), 
Schwartz D surface (middle), and Schoen G surface (right). 

Recently, there has been a resurgence of interest in triply periodic minimal sur- 
faces because two-phase composites whose phases are separated by such surfaces 
have been shown to be optimal with respect to multifunctional optimization of their 
material properties [El IS El- For instance, it has been shown that maximal 
values for the sum of the effective thermal and electrical conductivities is achieved 
by two-phase structures separated by Schwartz P and Schwartz D surfaces [ HH [T2] . 
These structures have also been discovered to be optimal for multifunctional bulk 
modulus and electrical conductivity optimizations [|l3j. As another example, the 
porous medium with the Schwartz P interface was found (via computer simulation) 
to have the largest fiuid permeability over a range of microstructures examined [ E] ■ 

In [ Jung and Torquato observed that the fiuid permeabilities of porous, 
bicontinuous microstructures is inversely related to the total interfacial surface area 
per unit volume. This observation led them to conjecture that the maximal fiuid 
permeability for a triply periodic porous medium is achieved by the structure that 
minimizes the total interfacial area One of the most interesting aspects of 

this conjecture is its focus on a global property of the interface {i.e., total surface 
area) rather than on local {e.g., differential) properties of the surface {e.g., mean 
curvature) . 

terms of the divergence of the unit normal vector (or equivalently, the trace of the gradient of n): 
H{-k) = -iV • n = -itr(gradn). 
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Unfortunately, exploration of global properties of surfaces, such as the total sur- 
face area, seems to have received little attention in the literature. Classically, the 
study of surfaces and their properties has been the realm of differential geometry. 
In that field, research has traditionally focused on characterization and examina- 
tion of minimal surfaces {i.e., surfaces with zero mean curvature). The search for 
minimal surfaces has been ongoing activity since the mathematician H. A. Schwartz 
published the first example of a minimal surface with full three-dimensional period- 
icity in 1865 [|15j. A systematic exploration of triply periodic surfaces with nonzero 
mean curvature appears to have only been carried out relatively recently by An- 
derson, Davis, Scriven and Nitsche [Ej. In their work, Anderson et al. numerically 
computed surfaces of prescribed mean curvature (possibly non-constant) and studied 
their properties. 

In this paper, we consider triply periodic surfaces of nonzero mean curvature from 
a slightly different, but complementary, perspective. Motivated by the fluid perme- 
ability conjecture and the scarcity of mathematical results on surfaces of nonzero 
mean curvature, we examine triply periodic surfaces that minimize the total interfa- 
cial surface area while satisfying a constraint on the volume fractions of the regions 
that the surface separates. Interestingly, we find that the optimal surfaces for this 
problem are precisely those possessing a constant mean curvature. This result al- 
lows us to reproduce many of the results obtained by Andersen et al.. However, 
we emphasize that our approach comes from a fundamentally different perspective. 
Whereas Andersen et al. prescribe the mean curvature and consider properties of the 
resulting surface, we specify the volume fraction and ask what the mean curvature 
of the surface must be to minimize the total surface area. In general terms, we are 
interested in the properties of structures that arise when global geometric proper- 
ties of a surface are prescribed (as opposed to examining the consequences of local 
properties). 

The main conceptual foundations for our work come from the fundamental ideas 
underlying the level set method [ HHI HZj . Using these ideas, we have developed a 
novel numerical method for computing surfaces for which the surface local 
minimum and that satisfy a specified volume fraction constraint. The ideas underly- 
ing the level method were also used to provide a theoretical characterization of the 
surfaces that solve of our constrained optimization problem. For both computational 
and theoretical purposes, the level set method formulation is convenient because it 
does not require explicit representation of the surface. From a theoretical perspec- 
tive, the level set method also has the advantage that the formulas for key geometric 
quantities, such as the surface area and volume, have a relatively simple form. 

In section |21 we begin by developing a level set formulation for describing and 
analyzing triply periodic surfaces. In section|2l we use this framework to theoretically 
characterize extrema of the total interfacial area under a constraint on the volume 
fraction of the phases. In sectional we present a numerical optimization procedure 
for computing surfaces that are local minima of the total interfacial area subject 
to a volume constraint. Finally, in section we use a parallel implementation of 
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our numerical algorithm to study several physically and mathematically interesting 
surfaces. Parallelism was utilized to help deal with the computational cost of our 
fully three-dimensional optimization problems. 

2. Level Set Formulation 

Our exploration of triply-periodic surfaces that are local minima of the total 
interfacial area subject to a volume constraint is based on ideas underlying the ap- 
plication of the level set method to shape optimization [UHllin]- Following the usual 
methodology E], we represent surfaces as the zero level set of an embedding 
function, (f){x), defined throughout the three-dimensional volume and use variational 
calculus to derive the relationship between (j) and global geometric properties of the 
surface. 

In this section, we present a detailed discussion of the level set formulation for 
surface area minimization in the presence of a volume fraction constraint. While 
it is clear that similar formulations have been used in previous studies [ |]^ , the 
explicit formulas and derivation for the most important geometric quantities were 
not provided. In the present discussion, which is intended to fill this apparent gap in 
the literature, we derive formulas for the area and volume of a triply-periodic surface 
and the variations of these quantities with respect to (p. It is worth mentioning that 
the following derivation also holds for non-periodic surfaces and is easily extended 
to codimension-one surfaces in n-dimensional Euclidean space. 

For the triply periodic problems that we are interested in, the domain is taken 
to be a unit cell of the periodically repeating structure. Let the surface of interest, 
r, be represented by the zero level set of (p. Then F divides the unit cell into two 
distinct phases. Without loss of generality, we define the region where 0(x) < to 
be phase 1 (see Figures El and . An important assumption in our work is that a 
shift of the unit cell by a lattice vector does not cause an interchange of the phases. 
In this paper, systems possessing this property will be called phase-periodic. 

In terms of the embedding function, 0, our minimization problem may be stated 

as 

minimize ^(0) subject to /(0) = fo, (1) 

where A{4>) is the total surface area of the zero level set, /(0) is the volume fraction 
of phase 1, and fo is the desired volume fraction for phase 1. Unfortunately, there 
are currently no theoretical results concerning the global minimum of our problem. 
However, as we shall see, there are many local minima for this problem. Because 
many interesting surfaces arise as local minima, we shall focus our attention on un- 
derstanding the structure of these surfaces and neglect, for the moment, the question 
of global minimality. 

We use the method of Lagrange multipliers to solve the optimization problem (^J. 
The Lagrangian is given by 



£(0, A) = ^(0) + A (/(0)-/o), 



(2) 
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Figure 2. Schematic diagram of a codimension-one surface and the sign of the level 
set function, (p, in various regions of the domain. The shaded and unshaded regions 
correspond to phase 1 and phase 2, respectively. 




Figure 3. The phase 1 region corresponding to three different choices for F: Schwartz 
P surface (left), Schwartz D surface (middle), and Schoen G surface (right). 



where A is the Lagrange multiplier. Taking the variation of C with respect to 0, we 
find that a necessary condition for a minimizer is 

6C{(P,\) = 6A{(f)) + \6f{(P) = 0. (3) 
This condition, together with the volume fraction constraint 

/(0) - /o = (4) 
allows us to compute (p and A. 

2.1. Volume integral formulation of surface integrals 

To compute the variations of A and / required by our Lagrangian formulation of 
the optimization problem, it is convenient to first derive a relation that allows us to 
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convert between surface integrals over T and volume integrals over entire unit cell. 
Let Q denote the entire unit cell, and let N denote an outward pointing unit normal 
vector on the boundary of the unit cell, dQ. Similarly, let fli and n denote the region 
corresponding to phase 1 and an outward pointing normal vector on F = dQ,i. Now, 
consider the surface integral of an arbitrary function over F: 



p{4>) dS. (5) 

Applying the divergence theorem, we can rewrite this surface integral as a volume 
integral over phase 1 and surface integrals over the intersections of phase 1 with the 
boundaries of the unit cell: 

J p{(j)) dS = J p(0)n -ndS = j ■ n dS 

where B is the union of the intersections of phase 1 with the boundaries of the 
unit cell and we have identified n with n^^- Notice that the last surface integral 
term vanishes because we have assumed that is phase-periodic which means that 

j is invariant under a translation of x by any lattice vector. Using phase- 
periodicity, we can cancel out contributions of surface integrals from opposite faces 
of the unit cell because N on opposite faces have opposite sign. 

Next, we observe that the volume integral over phase 1 can be rewritten as a 
volume integral over the entire unit cell: 

where 0(0) is the Heaviside function 



0(0) 



> 

< 0. 



This expression can be further simplified via an integration by parts procedure to 
obtain: 



p{<P)dS = -^V(l-0(0))-(^p(0)^)ciy 



+ / (1-0(0))(M0)^)-Nd5 



= / p(4>) S(<t>) IIV0II dV. (8) 
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where 5{(j)) is the Dirac delta function and we have ehminated the surface integral 
contribution by appealing to a similar line of reasoning as used when eliminating 
the surface integral term from The form of (jH)) is not unexpected - the delta 
function ensures that only values of p{(f)) on the surface F contribute to the integral 
and the ||V0|| is the Jacobian that arises when (as opposed to a spatial coordi- 
nate in the direction normal to the surface) is used as the argument of the delta 
function. Equation (jH} forms the foundation for much of the following mathematical 
development. 



Jj^dS, and (jHl), we find that 



2.2. Variation of the surface area 

Using the formula for the area of F, A{(f)) 

A{<j>) = [ 5(0) IIV0II dV 
Jn 

Taking the variation of this expression with respect to yields 

5^(0) = /[5{5(0)}||V0||+5(0)5||V0||]dl- 
Jn 

V0- V(50)" 



(9) 



5'(0)||V0||50 + 5(0) 

V0 



V(5(0)) 



50 + 5(0) 



liv^ll 

V0- V(50) 



dV 



(10) 



IIV0II " IIV0II 

where 5'(0) is the first derivative of the Dirac delta function and we have used the 
relationship V5(0) = 5'(0)V0. We can further simplify 6A{(f)) by using the product 
rule to make the following substitutions: 

V0 



V(5(0)) 



and 



5(0) 



IIV0II 

V0- V(50) 
IIV0II 



V- 



50V 



(11) 



(12) 



Thus, (II U|) can be rewritten as 



5^(0) 



V 



505(0)V ■ f 



VIIV0||J 



dV 



(13) 



where we have again used phase-periodicity to eliminate the surface integral over 
the unit cell that arises when applying the divergence theorem to the first integrand. 
Finally, we may convert 6 A back into a surface integral over F by comparing (jl3j] 
with (B: 



5^(0) = -^(V-n) 



50 



dS. 



(14) 
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2.3. Variation of the volume fraction constraint 

To compute the variation of the volume fraction constraint, we begin by writing 
the definition of the volume fraction as an integral over the entire unit cell: 



/(0)= / idv= {1-emdv. (15) 

Jni Jn 

Note that in writing (jl5|) . we have implicitly assumed that the volume of the unit 
cell is 1. Next, we compute the variation of f{(j)) as follows: 

5/(0) = [ 6{1- 0(0)) dV = - [ 6{<j>)6<j>dV = - ^ ^^^^ 



where we have again made use of (jH)) to convert the volume integral over the entire 
unit cell to a surface integral over F. 

3. Theoretical Analysis of Minimization Problem 

3.1. Characterization of local extrema 

Using the machinery developed in the previous section, we can characterize the 
local extrema of the optimization problem (Q). A straightforward application of the 
calculus of variations with constraints leads to the following 

Theorem. ^(0) is a local extremum subject to the constraint f{(f)) = fo if and only 
if the surface T has constant mean curvature. 

Proof (^) Suppose that F is a local extremum of the total surface area subject 
to the specified volume fraction constraint. Then it is a local extremum (without 
constraint) of the Lagrangian [ 20;: 

£(0, A) ^ ^(0) + A (/(0) - fo) , (17) 

where A is a scalar Lagrange multiplier. Using the expressions for the variation of the 
surface area and volume fraction derived in the previous section, the total variation 
of £(0, A) is given by 

5£(0,A) = 5^(0) + 5 [A (/(0) - /„)] 

= - ^ [(V ■ n) + A] ^rf^ + (/(0) - /,) a. (18) 

Since this expression must be of zero for any variation in and A, (fTHj) implies that 

V-n + A = onF (19) 

/(0)-/o = 0. (20) 

Equation (j^n|) is just the volume fraction constraint (which is expected and arises 
for all applications of the method of Lagrange multipliers). Equation (fT^ . however, 
leads to the conclusion that the mean curvature is constant over the entire surface 
and is equal to half the Lagrange multiplier. 



Level Set Approach for Surface Area Minimization of Triply Periodic Surfaces 9 



(<^=) Now, suppose that F is a constant mean curvature surface. Then, the vari- 
ation in the area simphfies to yield 

= - ^(V . „)^dS = -(V . n) ^ = (V . nWW- (21) 

Since any variation of the level set function that does not change the volume 
fraction must satisfy 5/(0) = 0, we see that 6A{(p) = for any volume fraction 
preserving variation in 0. In other words, F is a local extremum of the total surface 
area subject to a volume fraction constraint. □ 

Before moving on, it is worth mentioning some special cases of local extrema of the 
total surface area with a constrained volume fraction. First, there are the minimal 
surfaces. Because, they have zero mean curvature, the theorem immediately implies 
that all minimal surfaces are local extrema of the total surface area under a volume 
fraction constraint. However, minimal surfaces also have the interesting property 
that they are local extrema of the total surface area without any constraint. This 
result follows directly from the expression for the variation in the total surface area 
because the integrand is identically zero for minimal surfaces. A second class of 
important examples is the class of constant, nonzero mean curvature surfaces. Two 
examples are the sphere and infinite cylinder. Both of these objects have a constant, 
nonzero mean curvature, and thus are local extrema of the total surface area when 
the volume is constrained. Further examples are provided by any of the nonzero 
mean curvature surfaces studied by Anderson et al. [ITU]. 

These examples demonstrate that the class of minimal surfaces and the class 
of local extrema of the total surface area under a volume fraction constraint are 
truly distinct. That surfaces arising from an analysis of local geometric properties 
differ from those arising when global geometric properties are studied underscores 
the importance of considering the global geometric features of surfaces. These con- 
siderations are especially important from a physical perspective; while the driving 
force for surface evolution is often local in nature, global constraints {e.g., conserved 
quantities) are almost always present and can affect the global structures that arise. 

4. Numerical Optimization Procedure 

While Theorem 13. II provides a theoretical characterization of local minima of the 
total surface under a constrained volume fraction, it is does not provide a means 
for obtaining optimal surfaces. In this section, we present a numerical procedure 
for computing locally optimal surfaces. Our approach follows the work of Osher 
and Santosa [ 18j which evolves the embedding function along steepest descent 
directions using the evolution equation: 

50 + t;(x)||V0|| = O, (22) 

with f (x) chosen so that remains in the set of embedding functions that satisfy 
the volume fraction constraint. As pointed out in [d], this equation is equivalent 
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to a Hamilton- Jacobi equation if the change in is viewed as occurring continuously 
in time. In addition, we use an auxiliary Newton iteration to explicitly enforce 
the volume fraction constraint when numerical error causes the volume fraction to 
drift beyond an acceptable tolerance and to generate initial structures that satisfy a 
prescribed volume fraction constraint [E]- 

4.1. Descent direction and projected gradient algorithm 

Following [ilSj, we choose the velocity field f (x) to be the steepest descent direc- 
tion for the Lagrangian, C From ()18|1 . we know that the variation in C is given by 



5C{^, A) = - / [(V ■ n) + A] ^^dS. (23) 



r 



Therefore, choosing 50 such that 

50= [(V-n) +A]||V0|| (24) 



ensures that C decreases at every iteration. Comparing with (jzzj) . we can identify 
the velocity field w(x) as 

t;(x) = - [(V ■ n) + A] . (25) 

This velocity field has a natural interpretation: the motion of the interface is driven 
by mean curvature (which tends to shrink the local area of the interface) and the 
Lagrange multiplier (which tries to keep the volume fraction from changing). 

To compute the value of A to use for each iteration, we insist that 50 is chosen 
to satisfy 

5/(0) = 0. (26) 

Essentially, we are ensuring that + 50 satisfies the constraint equation linearized 
about the current iteration of 0: 

/(0 + 50)^/(0) + 5/(0) (27) 

Substituting into (fTB|) and enforcing we find that the Lagrange multiplier 
should be set equal to 

_ J,{V-n)dS _ /r(V ■ n)dS 
^~ J^ldS ~ Ai<P) ■ ^^^^ 



Notice that A is equal to twice the average value of the the mean curvature over the 
surface. 
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4.1.1. Extension of the velocity off of the interface 

The velocity of points off of the zero level set are obtained by extending the 
velocity on F in the normal direction [Ej. We decided to set the velocity in this 
way to avoid the mathematical singularities that arise for non-zero level sets with 
very large mean curvatures when using the "natural" velocity extension^. These 
singularities, which appear at cusps of the embedding function 0, lead to numerical 
difficulties [e.g., strict stability constraint on the effective time step size). 

4.2. Newton iteration to enforce volume fraction constraint 

Because we are only approximately enforcing the constraint at each time step, 
the iterates of the embedding function will eventually fail to satisfy the volume 
fraction constraint. To put an iterate back onto the feasible set after the constraint 
has been violated by more than an acceptable tolerance, we use a Newton iteration as 
in [CHI- ^^^o the Newton iteration to ensure that the initial condition for the 
embedding function satisfies the desired volume fraction constraint. For situations 
where the target volume fraction constraint is far from the actual volume fraction 
of the initial iterate, we use continuation in fo to improve the convergence of the 
volume fraction constraint algorithm [ 21]. 

We begin by considering the volume fraction of a corrected embedding function: 
+ 50), where (p^^^ is the uncorrected embedding function. Next, we think of 
S(f) in ()24j) as a function of A (taken to be an unknown) and introduce a scale factor 
a > 0: 



This choice for 50(A) was originally proposed by Osher and Santosa in [ 18 . While 
there are certainly other possible choices for 50(A), we opted to use because it 
tends to preserve the shape of the surface, it performed well in the present work, and 
it allowed us to reuse code written for other portions of the computation. 

Note that ^ allows us to think of /(0(°) + 50) as a function of A. We can then 
use a Newton iteration to compute the value of A (and therefore the appropriate 
correction 50) required so that the volume fraction constraint is satisfied: f{4>^^^ + 
50) = 0. It is important to recognize that the choice 50 in the Newton iteration is 
different from the choice of 50 in the optimization iteration (j23). In the latter case, 
A has a known value and can be explicitly computed; in the former case, A is an 
unknown variable that is determined through a Newton iteration. 

Taking advantage of the level set formulation of our problem, calculating the 



^The "natural" velocity extension uses H25|) for points off of the zero level set with n replaced by 

V0 




(29) 



iiv</.ir 
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derivative of /(A) with respect to A is straightforward: 

Da/(0^°) + 50(A)) = - f D^H{(P^'^^ + 6(P)dV 

Jn 

= - /" D^(O)+,^i/(0(°)+50)DA(0^°)+50)c?V^ 

Jn 

= -a [ (5(0(°) + 50) ||V0||(il^. (30) 
Jn 

4.2.1. Choice of scale factor 

The scale factor is necessary to ensure that the zero level set of (^(j)^^^ + 50) does 
not vanish and remains within the unit cell; otherwise, the derivative of / with respect 
A which is computed via as integral over the unit cell in (j3n|l may vanish. Intuitively, 
the scale factor helps make the discretized version of 50 a reasonable approximation 
to an infinitesimal variation of 0. Without it, the discrete approximation to 50 would 
be equal to [(V • n) + A], which may not be small for all Newton iterations. In our 
simulations, we employed a scaling factor of a = (min{Axi, Ax2, Ax^})"^ where Axt 
is the grid spacing in each i-th coordinate direction. This choice for a was selected 
by numerical experimentation. 

4.3. Numerical implementation issues 

The optimization algorithm based on the algorithmic components discussed in 
the previous two sections may be summarized in the following steps: 

1. Generate an initial configuration for 0o(x). 

2. Evolve the level set function to time tn+i if > tol. 

(a) If necessary, use Newton iteration to enforce the volume fraction con- 
straint. 

i. Select an initial guess for A. 

ii. Repeat the following steps until convergence: 

A. Compute 50 using (|^. 

B. Compute 0^°) + 50. 

C. Update A using: A^'^+i) = A^^) - ^j^^)('^^oT|g^(AW)^° ^^^''^ ^ 
A;-th Newton iterate. 

(b) Compute the Lagrange multiplier A„ using (j2Hl)- 

(c) Compute the descent direction 50„ using (j2D). 

(d) Update using 0„+i(x) to 0n(x)+/350n(x). /3 is a scale factor required for 
stability of the numerical scheme. It is equivalent to a stable time step size 
when (|22j) is viewed as a Hamilton- Jacobi equation varying continuously 
in time. 

(e) Periodically reinitialize 0(a;) to an approximate distance function within 
a sufficiently wide band around the zero level set of 0(x). 
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4.3.1. Synergy of reinitialization and volume fraction constraint enforce- 
ment 

Both reinitialization and enforcement of the volume fraction constraint introduce 
errors into (p. However, the two algorithms often work together in concert - each 
algorithm reducing the error introduced by the other. 

It is well-known that the reinitialization procedure can cause the zero level set of 
to shift from its "true" position [^lEj- Fortuitously, the process of enforcing the 
volume fraction constraint helps to reduce the magnitude of this shift by disallowing 
large shifts in the zero level set. On the other hand, the volume fraction constraint 
algorithm can cause a significant shift of the zero level set (especially when gener- 
ating initial conditions of a specified volume fraction). If a narrow banding type 
of procedure is used, the zero level set could end up in a region where is not 
very smooth. Reinitializing in a sufficiently wide band prior to enforcing the volume 
fraction constraint helps to avoid this problem. Reinitialization after enforcing the 
volume fraction constraint also helps to keep the zero level set "centered" within the 
narrow band that is actively updated by the optimization algorithm. 

4.3.2. Stopping criteria for optimization loop 

Because our goal is to minimize A, the change in the interfacial surface area 
between iterations, AAn is used as the stopping criteria. Thus, the embedding 
function, 0, is updated until the AAn falls below a prescribe tolerance. 

4.3.3. Evaluation of surface integrals 

We numerically evaluate surface integrals by first converting them to volume 
integrals via (jHI) and then replacing the delta functions with smoothed approximates. 
For example, using (jH)), the Lagrange multiplier (j^H|) may be written in the form 

J^{V-n)6{<P)\\V<P\\dV 

J^6i<l^)\\V<l>\\dV ■ ^ ^ 

In our computations, we use the following approximation for the 5- function [I16j: 




^^^^^ = , ,.ii<^iKi r;; J ' (32) 

cos(^)] II0II < e, 

where e is the width of the smoothed 5-function. Following the usual conventions for 
level set method calculations, we choose e = 3 max{Axi, Ax2, Axs}. 

4.3.4. Discretization of the mean curvature 

The discretization of the curvature term is performed by using the natural gen- 
eralization of the first-order accurate scheme described by Zhao, Chan, Merriman, 
and Osher 



4.3.5. Grid resolution requirements 

All of our numerical solutions were computed on a unit cube using uniform meshes 
with sizes between 100 x 100 x 100 and 250 x 250 x 250. For relatively simple surfaces, 
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such as the Schwartz P surface, the lowest resolution mesh was sufficient to obtain 
accurate values for the total surface area^. However, complex surfaces, where the 
distance between separate sheets of the surface is smaller, required higher resolution 
meshes. For example, in order to obtain results comparable in accuracy to those 
for the Schwartz P and Schoen G surfaces with a mesh of size 200 x 200 x 200, the 
Schwartz D surface calculations required a grid of size 250 x 250 x 250. Finer meshes 
were also required to obtain accurate values for the mean curvature and when finding 
optimal surfaces for volume fractions other than 0.5. 

4.3.6. Parallel computation 

Due to the high computational cost (both in time and memory) of the fully 
three-dimensional shape-optimization calculations, we implemented our algorithm 
using LSMLIB [^5, a parallel level set method software library developed by one of 
the authors. Rapid implementation of our parallel simulation code and LSMLIB was 
achieved by leveraging the parallel computing framework built into the Structured 
Adaptive Mesh Refinement Application Infrastructure (SAMRAI) [EHESj developed 
at Lawrence Livermore National Laboratory (LLNL). 

5. Numerical Results 

In this section, we present results obtained using the numerical optimization 
algorithm developed in Section |3] Simulations were run on Linux clusters^ using 
between 2 and 8 nodes depending on the size of the computation. Results were 
visualized using the Visit program developed at LLNL. 

5.1. Verification of numerical optimization algorithm 

To verify our numerical optimization algorithm, we tested its ability to accurately 
compute a few structures whose local optimality are easily verified. These results 
were computed on a grid of size 100 x 100 x 100. 

Our ffist test case is the cylinder. Because it has constant mean curvature. 
Theorem 13.11 shows that it is a local minimum of the total surface area under a volume 
fraction constraint. From symmetry considerations, we expect that any infinitely 
long channel should evolve towards the locally optimal cylinder structure. In our 
numerical calculations, we started the optimization algorithm with an infinitely long 
square channel as the initial surface. As expected, the square channel evolves to the 
optimal cylindrical channel (see Figure |3)). 

Our second test case is the sphere. Again, Theorem 13.11 shows it is a locally 
optimal surface. In this test, we started the optimization algorithm with an closed 
cube (chosen from symmetry considerations) as the initial surface. Figure El shows 
that the cube evolves to the sphere, as expected. 

■^The difference in the total surface area computed at the 100 x 100 x 100 and the 150 x 150 x 150 
resolutions was less than 10~*. 

■*Runs were performed on either a 75 node dual Opteron cluster or a 128 node Intel P4 system. 
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Figure 4. Evolution of an infinitely long square channel (left) to the locally optimal 
cylindrical channel (right). 




Figure 5. Evolution of a cube surface (left) to the locally optimal sphere surface 
(right). 

In both calculations, the volume fraction constraint algorithm was used to enforce 
fo = 0.5. However, the volume fraction did not fluctuate much during the numerical 
simulations, so the volume fraction constraint algorithm was rarely invoked. Also, 
notice that for both of these test cases, the level set calculations do not appear to 
have any difficulty handling sharp edges/ corners. 

5.2. Local minimality of Schwartz P, Schwartz D, and Schoen G surface 
areas 

There are many known ways to compute the Schwartz P, Schwartz D, and Schoen 
G triply periodic minimal surfaces. They can be characterized exactly using an 
Enneper-Weierstrass (complex integration) representation [ generated as the 
local minima of the scalar order parameter Landau- Ginzburg functional used to 
describe ordering phenomena in microemulsions [I2Z|; and approximated by Fourier 
series using the periodic nodal surface (PNS) expansion [ |2H1 12 123 EO] • It is worth 
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pointing out that the approach of minimizing the Landau-Ginzburg functional is 
actually a phase-field version of our approach for an appropriately chosen form of 
the energy functional. 

Because minimal surfaces have zero mean curvature everywhere, Theorem 13.11 
guarantees that the Schwartz P, Schwartz D, and Schoen G surfaces are locally opti- 
mal surfaces. However, it does not indicate whether these surfaces are local maxima 
or minima or saddle points of the total surface area. The numerical optimization 
procedure developed in Section 14.31 provides a means to answer this question. Since 
the optimization procedure seeks surfaces with minimal total surface area, locally 
minimal surfaces should be stable under perturbations. In contrast, local maxima or 
saddle points of the total surface area should be unstable^. 




Figure 6. Locally optimal surfaces for /o = 1/2 starting from the PNS approximation 
to the Schwartz P surface (left), Schwartz D surface (middle), Schoen G surface 
(right). Notice that the final surfaces are precisely the Schwartz P, Schwartz D, and 
Schoen G surfaces, which indicates that all three of these surfaces are local minima of 
the total surface area when the volume fraction is 1/2. These results were computed 
using a grid of size 150 x 150 x 150. 



In our numerical simulations, we start the optimization procedure using the PNS 
approximations^ (which are very good) as the initial (j). As can be seen in Figure 
the optimization procedure evolves to the correct triply periodic minimal surface. 
Therefore, we can conclude that the Schwartz P, Schwartz D, and Schoen G surfaces 
are all local minima of total surface area when fo = 1/2. 

^Note that instability is a necessary by not sufficient condition for a local maximum or saddle point 
while stability is a sufficient by not necessary condition for a local minimum. Thus, we can only 
safely draw conclusions when the surface is a local minimum of the total surface area. 
^In the PNS approximation, each term consists of sine and cosine functions where the high-order 
terms have a high reciprocal lattice vector norm [I29|. 
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5.2.1. Width of basins of attraction for Schwartz P, Schwartz D and 
Schoen G 

By starting the optimization procedure with other initial configurations of care- 
fully controlled symmetry, we were able to qualitatively explore the width of the 
basin of attraction for the Schwartz P, Schwartz D, and Schoen G surfaces. 

For the Schwartz P surface, we considered two different initial configurations: (1) 
triply periodic circular channels (see Figure [7j) and (2) set to a function dominated 
by the second term in the PNS approximation of the Schwartz P surface. As seen in 
Figure the triply periodic circular channels configurations evolves to the Schwartz 
P surface. However, the second initial configuration leads to a new high genus^ 
surface with the same symmetry as the Schwartz P surface (see Figure |H1). The 
symmetry of the computed locally optimal solution is expected from symmetry of 
the initial configuration, but the genus of the solution is a bit unexpected. Clearly, 
this result is related to the width of the basin of attraction for the Schwartz P surface. 



Figure 7. Evolution of triply periodic circular channels (left) to the Schwartz P 
surface (right). 



For the Schwartz D and Schoen G surfaces, we generated initial configurations by 
setting the initial condition for to be a function dominated by the second term in 
the PNS approximation of the Schwartz D and Schoen G surfaces, respectively. This 
procedure allowed us to preserve the symmetry of the of respective surface, while sig- 
nificantly perturbing the actual surface. Figure IHl shows that the Schwartz D surface 
is the locally optimal structure even if the second term in the PNS approximation 
is dominant. Even though the initial configuration in Figure El looks completely dif- 
ferent from the Schwartz D surface, it converges to Schwartz D surface. This result 
implies that the Schwartz D has a relatively wide basin of attraction. The Schoen G 



^The genus describes how many holes are in a closed surface and therefore is an integer number. In 
the infinitely periodic case, we describe the genus in a unit cell, which will be a finite number [ll(J|. 
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Figure 8. Evolution of a triply periodic surface with the symmetry of the Schwartz 
P surface to a new high genus minimal surface (right). The initial structure (left) 
is obtained by modifying the PNS approximation to the Schwartz P surface so that 
the second term dominates. 

surface shows a similar wide basin of attraction (see Figure irU]) . It should be noted 
that the three locally optimal structures all have zero mean curvature at every point 
on their surfaces. 




Figure 9. Evolution of a triply periodic surface with the symmetry of the Schwartz 
D surface to a the Schwartz D surface (right). The initial structure (left) is obtained 
by modifying the PNS approximation to the Schwartz D surface so that the second 
term dominates. 



5.3. Locally optimal structures when fo ^ 1/2. 

In this section, we examine three families of surfaces with the symmetry and 
topology as the Schwartz P, Schwartz D, and Schoen G surfaces but with different 
values for the volume fraction of phase 1. We generated initial configurations with 
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Figure 10. Evolution of a triply periodic surface with the symmetry of the Schoen G 
surface to a the Schoen G surface (right). The initial structure (left) is obtained by 
modifying the PNS approximation to the Schoen G surface so that the second term 
dominates. 



the desired symmetry, topology, and volume fractions by taking the PNS Schwartz 
P, Schwartz D and Schoen G approximations, respectively, and using the volume 
constraint algorithm to impose the desired volume constraint. We note that to obtain 
initial conditions for some of the volume fractions reported, we used continuation in 
the volume fraction to improve the convergence properties of the Newton iteration 
scheme. These initial configurations were then allowed to evolve to the nearest 
optimal structure. The total surface area A is computed using © and the mean 
curvature H is obtained from (jHlj) using the relation H = 0.5A. In both calculations, 
we have used the smoothed delta function, 6^, defined in (IH^ . 

5.3.1. Schwartz P surface family 

We carried out the shape optimization procedure for the following volume frac- 
tions: 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75. Figure CU shows the 
six different Schwartz P-type surfaces having specified volume fractions. These cal- 
culations were performed on a grid of size 200 x 200 x 200. For a typical simulation, 
the area stopping criterion AAn < 10~^ was reached after approximately 29 hours 
for a 4 node parallel calculations on the Opteron cluster. The computation time 
in general increases for volume fractions further from fo = 0.5 because the volume 
fraction constraint algorithm is invoked more often. 

As the Theorem I3.1l indicates. locally optimal surfaces should have constant mean 
curvature. Because each surface in the Schwartz P family is locally optimal, they all 
constant mean curvature. In Table Q we report the numerical values of the mean 
curvature and total surface area for the surfaces in the Schwartz P surface family 
that we computed. Figure IT^ compares our results with those obtained by Anderson 
et al. ^. Note that the data we produce matches the results of Anderson et al. rea- 




The curves for the data obtained by Anderson et al. are reproduced with the permission of the 
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Figure 11. Family of surfaces possessing the symmetry and topology of the Schwartz 
P surface but with different volume fractions for phase 1. Upper left to upper right 
panel: fo = 0.25,0.35,0.45. Lower left to lower right panel: fo = 0.55,0.65,0.75. 

Table 1 

The mean curvature and total surface area of the Schwartz P surface family for 



different values of the volume fraction. 


Volume fraction fo 


H 


A 


0.25 


1.67 


2.00 


0.3 


1.12 


2.14 


0.35 


0.77 


2.23 


0.4 


0.50 


2.30 


0.45 


0.24 


2.33 


0.5 


0.00 


2.34 


0.55 


-0.24 


2.33 


0.6 


-0.49 


2.30 


0.65 


-0.78 


2.23 


0.7 


-1.12 


2.14 


0.75 


-1.67 


2.00 



sonably well. The small discrepancy between the two sets of data is a consequence 
authors and were generated by digitizing the plots in [JMi using the Plot Digitizer program. 
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of (1) the lower effective grid resolution of the calculations by Anderson et al. for 
the Schwartz P surface and (2) the digitization procedure used to extract the curves 
of Anderson et al. from a printed version of their paper. Figure ^1 shows that the 
volume is almost linearly related to mean curvature in the volume fraction range be- 
tween fo = 0.35 and fo = 0.65. Beyond this range, the mean curvature varies more 
rapidly as a function of the volume fraction. Anderson's data shows that near volume 
fractions of 0.75 and 0.25, the relationship between the volume fraction and mean 
curvature reverses direction. Due to the limitations of our computational method, 
we have not yet been able to reproduce the portions of curve between points B and 
C. We conjecture that our numerical method fails for these portions of the curves 
because they are not local minima of the surface area. 




Mean Curvature Volume Fraction 



Figure 12. Relationships between the volume fraction, mean curvature, and total 
surface area for optimal surfaces of the Schwartz P family: volume fraction versus 
mean curvature (left) and total area per unit cell versus volume fraction (right). Our 
results (circles) are compared with those obtained by Anderson et al. (solid lines). 



5.3.2. Schwartz D surface family 

The Schwartz D surface has largest surface area among the three triply peri- 
odic surfaces studied in this paper. Since large surface area is inversely correlated 
with minimum distance between distinct sheets of a surface, the calculations for the 
Schwartz D surface required higher grid resolution than those for the Schwartz P 
surface. For the present calculations, we employed a grid of size 250 x 250 x 250. 
Coarser meshes [e.g., 150 x 150 x 150) failed at extreme values of the volume frac- 
tions {e.g., fo = 0.15,0.85) because the optimal surfaces at these volume fractions 
have very narrow connections. Due to the larger computational grid, the run time 
for these simulations took significantly longer than for the Schwartz P-type surfaces. 
Reaching an area stopping criterion of AAn < 10~^ took approximately 30 hours for 
a 4 node parallel calculation on the Opteron cluster cluster. As with the Schwartz 
P-type simulations, the running time depended on the volume fraction constraints. 
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Figure 13. Family of surfaces possessing the symmetry and topology of the Schwartz 
D surface but with different volume fractions for phase 1. Upper left to upper right 
panel: fo = 0.15, 0.3, 0.45. Lower left to lower right panel: fo = 0.55, 0.7, 0.85. 

Figure IT^ shows excellent agreement with the results of Anderson et al. ^. How- 
ever, as with the Schwartz P-type surfaces, our algorithm cannot compute surfaces 
between points B and C on the curve. The numerical values of the mean curvature 
and total surface area for the surfaces in the Schwartz D surface family are reported 
in Table H 

5.3.3. Schoen G surface family 

For these calculations, we employed a grid of size 200 x 200 x 200. Because 
Anderson et al. did not study Schoen G surfaces, we explored the valid range of 
volume fractions by stepping the volume fraction in increments of 0.05 from a vol- 
ume fraction of 0.5 until we observed optimal surfaces that no longer possessed the 
symmetry of the Schoen G family fFigure IT3j). We detected the symmetry change 
(as well as sudden changes in the mean curvature and total surface area) beyond 
fo = 0.8 and / = 0.2. These results were consistent with simulations on finer meshes 
{e.g., 250 x 250 x 250), which leads us to conclude that the fo = 0.8 and fo = 0.2 
are close to the turning points in the volume fraction versus mean curvature graph 
(Figure HH)). 

As shown in Figure ^1 the mean curvature range for the Schoen G family of 



^In the calculations of Anderson et al., the effective grid resolution for the Schwartz D surface was 
higher than for the Schwartz P surface. 
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Figure 14. Relationships between the volume fraction, mean curvature, and total 
surface area for optimal surfaces of the Schwartz D family: volume fraction versus 
mean curvature (left) and total area per unit cell versus volume fraction (right). Our 
results (circles) are compared with those obtained by Anderson et al. (solid lines). 



Table 2 

The mean curvature and total surface area of the Schwartz D surface family for 



different values of the volume fraction. 


Volume fraction fo 


H 


A 


0.15 


3.86 


2.78 


0.2 


2.80 


3.11 


0.25 


2.11 


3.36 


0.3 


1.58 


3.54 


0.35 


1.12 


3.67 


0.4 


0.73 


3.76 


0.45 


0.35 


3.82 


0.5 


0.00 


3.84 


0.55 


-0.35 


3.82 


0.6 


-0.72 


3.76 


0.65 


-1.13 


3.67 


0.7 


-1.57 


3.54 


0.75 


-2.10 


3.36 


0.8 


-2.80 


3.11 


0.85 


-3.86 


2.78 



surfaces is similar to that of the Schwartz P family of surfaces in Section 15.3.11 
However, the total surface area of Schoen G-type surfaces are bigger than those of 
the Schwartz P-type surfaces with the the same volume fractions. Numerical values 
for the mean curvature and total surface area for surfaces in the Schoen G family 
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Figure 15. Family of surfaces possessing the symmetry and topology of the Schoen 
G surface but with different volume fractions for phase 1. Upper left to upper right 
panel: = 0.2, 0.3, 0.4. Lower left to lower right panel: = 0.6, 0.7, 0.8. 



are reported in Table El 



3.2 
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Mean Curvature Volume Fraction 



Figure 16. Relationships between the volume fraction, mean curvature, and total 
surface area for optimal surfaces of the Schoen G family: volume fraction versus 
mean curvature (left) and total area per unit cell versus volume fraction (right). 
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Table 3 

The mean curvature and total surface area of the Schoen G surface family for different 



values of the volume fraction. 


Volume fraction fo 


H 


A 


0.2 


2.15 


2.53 


0.25 


1.64 


2.72 


0.3 


1.22 


2.86 


0.35 


0.87 


2.97 


0.4 


0.56 


3.04 


0.45 


0.28 


3.08 


0.5 


0.00 


3.10 


0.55 


-0.28 


3.08 


0.6 


-0.56 


3.04 


0.65 


-0.87 


2.97 


0.7 


-1.22 


2.86 


0.75 


-1.64 


2.72 


0.8 


-2.15 


2.53 



Following Anderson et al. , it would be interesting to complete the graphs of the 
volume fraction versus mean curvature and total surface area. However, as with the 
Schwartz P and D surfaces, we are currently unable to generate surfaces beyond the 
turning points in the graph of the volume fraction versus mean curvature. 

6. Conclusions 

In this paper, we have generated and probed the structure of the constant mean 
curvature, triply periodic surfaces (including minimal surfaces) that arise when the 
total surface area is minimized subject to a constraint on the volume fraction of the 
regions that the surface separates. Specifically, we have shown that the Schwartz P, 
Schwartz D, and Schoen G minimal surfaces are local minima of the total surface 
area under a volume fraction constraint, studied the properties of families of surfaces 
possessing the symmetry of the Schwartz P, Schwartz D, and Schoen G surfaces 
at varying volume fractions of the phases separated by the surface, and generated 
new minimal surfaces (e.g. with high genus) as well as new surfaces with constant, 
nonzero mean curvature (e.g. family of surfaces with the symmetry of the Schoen 
G surface). Unlike many studies on surfaces, our perspective is interesting because 
it draws attention to the properties of surfaces when global geometric constraints 
rather than local differential conditions are imposed. Many of our results complement 
and extend the work by Andersen et al. on nonzero, constant mean curvature 
surfaces. We emphasize that an important feature of our approach is our ability to 
easily control the volume fraction of the phases. 

The basic approach that we have taken to study optimal triply periodic surfaces 
has been to use ideas from the level set method and shape optimization. This well- 
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known framework has allowed us to formulate the problem in a manner that is 
amenable both to theoretical analysis and numerical simulation. Theoretically, it 
allowed us to prove that surfaces with constant mean curvature exactly are precisely 
the ones that optimize total surface area under a volume fraction constraint on 
the phases. Computationally, the framework naturally led (via a steepest descent 
procedure) to a generalization of the two-dimensional shape optimization algorithms 
presented in [|^ and [EI- We note that our discussion of the optimization algorithm 
provides a more comprehensive description of how variational calculus is used in the 
context of shape optimization and level set methods. In particular, we demonstrate 
how the variational formulation of the problem can be used to draw theoretical 
conclusions (in addition to forming the foundation for the computational method). 

Using our numerical algorithm, we explored several optimal triply periodic sur- 
faces. When the volume fractions of the two phases are equal, our theorem implies 
that the Schwartz P, Schwartz D, and Schoen G surfaces are all local optima of the 
surface area; numerical simulation shows that these surfaces are actually local min- 
ima of the total surface area. For unequal volume fractions of the phases, we partially 
reproduced the volume fraction versus mean curvature and total surface area versus 
volume fraction results obtained by Anderson et al. [ QUI for the Schwartz P and 
D surfaces. Unfortunately, due to limitations in our algorithm, we were unable to 
reproduce the results for extreme values of the mean curvature. We conjecture the 
surfaces at extreme values of the mean curvature may no longer be local minima of 
the total surface area. This question certainly deserves require further investigation. 
In addition to studying the Schwartz P and Schwartz D families of surfaces, we also 
extended the work of Anderson et al. by determining the relationship between the 
volume fraction, mean curvature, and surface area for the Schoen G surface. 

One of the original goals our investigation was to show that the triply periodic 
surface with minimal total surface area at a volume fraction of 1/2 is the Schwartz 
P surface. While our numerical simulations provided us with empirical evidence 
supporting this conjecture, we could not prove these results because our numerical 
algorithm and theoretical analysis are limited to information about local minima of 
the total surface area. Thus, the question of the global optimality of the Schwartz P 
surface remains an open question for future investigation. 

6.1. Future directions 

To close, we mention a few future directions for our work. First, there are re- 
search questions that might allow us to gain more insight into the question of global 
minimality of the Schwartz P surface. For instance, a way to quantify the width of 
the basin of attractions for the local minima of the total surface area could allow 
us to be more systematic in our search through the space of triply periodic surfaces 
(which could be carried out using our numerical optimization algorithm.) Our nu- 
merical algorithm could also be used to explore the optimal surfaces possessing a 
wider array of symmetries. 

Second, with a working implementation of our numerical algorithm, we are now 
in a position to study important questions about the physical properties of two-phase 
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materials whose phases are separated by optimal surfaces. For instance, it has been 
conjectured that microstructures where two phases of unequal volume fraction are 
separated by a surface that minimize the total surface area lie on the upper bound 
of effective transport properties for composite materials [ El E]- One obstacle 
to collecting evidence supporting or refuting this conjecture has been our ability 
to generate surfaces of minimal surface area. Our optimization algorithm effectively 
removes this barrier and makes the multifunctional optimization problem for unequal 
volume fractions more tractable. 
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